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The Gilbert equation expressed in conventional polar angle 
notation is solved for one-dimensional Neel and Bloch walls in a 
ihin Permalloy film using different numerical methods to compare 
the extent of stability. The method.s examined include the Euler. 
so-called modified Dufort-Frankel and backward Euler methods. 
The backward Euler method is found to be stable under arbitrary 
magnitude of time difference for the case of Nee) walls if the contri- 
bution of the magnetization of the nearest-neighboring cells to the 
demagnetizing field was treated implicitly. The method is also 
stable in the case of Bloch walls, though the stability limit for the 
time step was only 10 times as large as that of the Euler method. 
The Dufort-Frankel method is found quite unstable. Two- 
dimensional compulation is found feasible using the backward Euler 
method. The principle and the results of the calculations arc given. 

Introduction 

While computation of wall structures in Permalloy films was 
mostly carried out based on variational method [ 1 1-(4], direct solu- 
tion of the Gilbert equation was reponed in abundance for domain 
vraXh in bubble materials. Because of the large magnetization of 
Permalloy, however, the equation includes (frastically increased 
demagnetizing torques, which may affect the stability of the numeri- 
cal solution. It is the purpose of the present papier lo examine 
whether or not the methods of discretizing the Gilbert equation 
which are useful in the case of bubble material are also useful in the 
case of Permalloy. 

Onc-dimcnsional Neel and Bloch wall structures arc calculated 
using the Euler [6], so-called modified Dufon-Frankel [51,(6 J and 
backward Euler [6] methods for that purpose. To examine the feasi- 
bility of large-scale computation , planar two-dimension wall struc- 
tures are also calculated using the backward Euler method. 

Definition of the problem 

We solve the Gilbert equation expressed in the polar and 
azimuthal angles, 6 and ^, of the magnetization vector Af as in [6] ' 
and [7J: 

(MAY)(ae -<) sine) = Sc/SO =/[0,(tij, 
(MiyXQ i- sin6) ^ (l/sine)5e/5(> =^ g 16,0), 

denoting with e the magnetic energy density of the system. 

The polar (z ) axis must be chosen in such a way that G never 
approaches zero (mod n) as described in [6] and [7\. Thus different 
coordinates systems are used to calculate Neel and Bloch structures 
(Fig. 1). While the main frame to be used for the Neel walls has the 
z -axis in the film normal, the Z-axis is in the film plane and perpen- 
dicular to the wall tangent in the subsidiary frame to be used for the 
Bloch walls. Using the main frame, for example, the magnetic 
energy density can be expressed as 

e[e. 0] = >l ((V8)2 + sin2e(70)2 1 - K sin^Osin^i^ -MHa 

denoting the self-demagnetizing field by H , 

In the one-dimensional calculations we consider a Pemnalloy 
film of thickness h and width xw (or Zw) as shown in Fig. 1. We 
divide the film along the a (or Z) axis into (or A'z ) prisms of 
equal size which extend in the y (orY) direction. The free bound- 
ary condition is used ai the film edges. 

In the two-dimensional calculations a rectangular film is 
divided into Nx by Ny cells of equal size as shown in Fig. 2. The 
easy axis is taken parallel toy - or 7 -axis. Using an initially defined 
angle parameter 5, the coordinates frame is selected at each cell (/, 
J) as follows: if the polar angle measured from the film normal is 
between n - 5 and 5 at cells (/, 7) and (/±1, /±1), the main frame 



is selected. Otherwise the subsidiary frame is selected. In Fig. 2 the 
subsidiary frame is used around the film center because of the 
existence of a Bloch line of "circulation" type. The parameter S is 
chosen between 10 and 80 degrees, depending on the mesh size 
used. The free boundary condition is used at all film edges. 

Because the purpose of the present paper is to compare the 
capability of different numerical methods, we make the following 
simplifications: 1) The Gilbert damping parameter is chosen arbi- 
trary. 2) Possible eddy current effects are not considered. 3) The 
magnetization is not allowed to vary through the thickness of the 
film. 

Numerical methods 

The numerical methods used to examine the one-dimensional 
calculation can be summarized as follows. In the Euler method the 
time derivatives are replaced by the respective forward differences. 
In the finite difference form of the reversible torques / and g the 
first and the second-order spaiial derivatives arc evaluated at the 
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Fig. 1 Coonlinalcs systems for (a) NccI anJ (b) Bloch walls. 
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Rg.2 Computing region and dual coordinates system for iwo- dimen- 
sional calculolion. An equilibrium state derived from h~ 300 A. 
&c s 5y e 70A, ATx* 8 and ^ = 12 is shown. 
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current time and discretized according to the midpoint fonnula. In 
ihc Dufort-Frankcl method the time derivatives arc discretized 
according to the midpoint formula. As to the spatial derivatives, a 
correction tcnn is added to the midpoint difference evaluated at the 
cuncnt time. In the present paper the method modified according to 
[51 so as to increase the stability is used [6]. 

In the backward Euler method the time derivatives arc replaced 
by the respective backward diflfcrcnce quotients with respective to 
time. The reference time is shifted to / 4-Af, one step after the 
current time, t . The lonjue tcmis / and ^ are also evaluated at 
/ + A/ . and expanded at time t so that temis of 0(A/X A9 and A^ are 
retained. In the Newton's method the Euler equation associated to E 
is solved, instead of the Gilbert equation, to obtain equilibrium 
configuration. The Euler ctjuation discretized according to the 
Newton's iteration scheme is nothing but the Gilben equation 
discretized according to the backwaid Euler method with time step 
A/ of infinity [6]. 

The spatial derivatives of exchange origin arc replaced by the 
linear combinations of the values at the nwsh points / under con- 
sideration and the two neighboring points /±l. Thus wc have a set 
of simultaneous linear equations in which AO and A^ at mesh 
points / and / ±1 are coupled t6).l7]. 

The demagnetizing field H is calculated with a finite grid 
method [IJ collecting all charges appearing on the surfaces of the 
cells within the computing region. Theoretically the backward Euler 
and the Newton's iteration schemes require to collect the contribu- 
tions of AO and A0 among all the mesh points in the computing 
region, which makes all the cells in the computing region coupled in 
the aforementioned simultaneous equarions. As in 161 and [7], how- 
ever» only the contribution of the cell under consideration and the 
neighboring cells is included in most of the present calculations so 
that we may not have simultaneous equations with a "dense" 
coefficient nutrix. Thus it would be more proper to call the 
presently used methods "incomplete Newton's method" and ''incom- 



<P (degree) 



h = 30, 100, 300 A 




plete backwaid Euler method" because long-range magnctostatic 
couplings are neglected. 

The actual calculation reveals that 1) Unlike the case with bub- 
ble material the Oufort-Frankel method is quite unsuible: 2) The 
incomplete backward Euler method allows us to use a time step at 
least 10 times larger than the largest time step which can be used 
with the Euler method; and that 3) In the case of Necl walls 
Newton's iteration schemes, both original and incomplete, converge 
stably, which would mean that the backward Euler method is stable 
with arbitrary magnitude of Ar . 

In the light of the results of one-dimension calculation the 
two-dimension calculation was performed using the incomplete 
backward Euler method only. The numerical techniques for two- 
dimension calculation such as the use of dual coordinates frame and 
implicit treatment of the demagnetizing field are the same as in [6] 
and [7] except for the use of the so-called ILUCXjS method [8] in 
the solution of die simultaneous linear equations. 

Results 

Wc use as material parameters M=800G. A =1.0x10"^ erg/cm, 
A:=1000erg/cm\ l7l=1.76xlO'iad/Oe-s and a^l.O. The equili- 
brium states obtained from the one-dimensional calculation are 
shown in Hg. 3. The initial value for 0 was chosen as 

= + 0.3jc(/ /Nx ). I ^Nxa 

= n/2 - 0^n{Nx-!)lNx . I>Nxa 

in the case of Necl walls. A step-like O jumping from -tU2 to 7i/2 at 
the wall center was used in the case of Bloch walls. A constant 
polar angle of lUl was used in both Neel and Bloch walls. The 
incomplete Newton*s niethod was used in the case of Neel walls. 
The calculation was repeated about 150 times until the rms value of 
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Fig.4 Average residual torque vs. calculation step derived from different 
methods. The case of Noel wall in a lOOA-thick film is shown. 
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Fig.3 Neel (a) and Bloch (b) walls in equilibrium. 
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Fig.5 Comparison of numencal stability. The configuration derived fmm 
140 iicra lions in the Newton's method (Fig.4) is used as the cora- 
inon initial state. 



3113 




35 

Oufort-Frankel 



10-6 U.4 X/t 

At«1.93x1(r^ns step 80 

F1g.6 Example of unstable calculation. 



the residual reversible torque decreased to 10~* of the initial value. 
The Bloch structures were calculated using the incomplete backward 
Euler method with At - 0.02ns. The calculations were observed to 
converge after - 100 steps. 

The speed of convergence observed in the calculation of the 
Neel wall using different methods is summarized in Fig. 4, where 
the used is indicated. The figure shows that incomplete 
Newton's iteration converges much faster than the Euler. or forward 
difference, method. We know, however, that incomplete Newton's 
iteration converges almost in several steps in the case of one- 
dimensional charged and Bloch walls in bubble matciial. The 
difference in the convergence speed would be ascribed to the 
extremely large demagnetizing field in Permalloy films. Although 
not practical, it is interesting to note that the original Newton's 
method converges in a few steps even in the case of Permalloy. For 
one-dimension Bloch walls in Permalloy films, which are accom- 
panied by much larger demagnetizing field, incomplete Newton's 
method was observed to diverge and the upper bound of ^ in the 
backward Euler method was only 10 times as large as that in the 
Euler method. The calculation using the original Newton's method 
was observed to converge even in that case. 

The modified Dufort-Frankel method is attractive in that it 
leads to an explicit solution not requiring solving simultaneous 
equations. The method was more stable than the Euler method in 
this case of bubble material [6]. The method, however, was found 
quite unstable when applied to Pfcrmalloy films as shown in Figs. S 



and 6. When M is varied, the method was stable if Af < lOOG. Fig- 
ure 6 shows a transient Necl type configuration which appeared in 
the process of an unstable solution obtained with M- lOOG. 

Figure 7 shows an example of planar two-dimensional calcula- 
tion performed with the backv^ard Euler method. The inplanc com- 
ponent of M was given initially in a concentric configuration around 
the film center. Ttit initial value of the film-normal component of 
M was given by 

6(7. y) = 2 arctan[cxp(4r(/, /)//b)] - tU2 

with and r denoting K{An%M^f^ and the distance of the mesh 
point (/ , J ) from the film center, respectively. The figure shows the 
configuration after 133 step calculation, where the rms value of the 
residual torque decreased to 5.4x10"' of the initial value. The com- 
pulation time was 10.75 minutcs/sicp on a MELCOM-COSMO 800 
with an operating speed of 1 .2 million instructions/s. The required 
main memory space was 1.6 Mbytes. 

Conclusion 

One-dimensional wall structures in thin Permalloy films were 
calculated using the Gilbert equation discretized according to vari- 
ous numerical methods in order to compare the capability of each 
method to suppress the numerical instability. The results of the cal- 
culation show that the backward Euler method gives the largest time 
step which can be used without inducing instability and that the so- 
called modified Dufon-Frankel method cannot be used in die calcu- 
lations for Permalloy filnis. The backward Euler method combined 
with ILUCXiS method is also shown to be useful for medium to 
large scale two-dimension calculation. On closing the authors 
would like to thank Prof. Takao Suzuki of Tohoku University 
(presently with IBM Research Division, Almaden Research Center) 
for his valuable comments on the present work. 
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Fig.7 Small closuic-Oux equilibrium conGgunitiQn derived from 2D planar calculation. 

The parameters used ait: A«20QA. 6xs5y»50A, Nx>30. Ny=60, and&fs0.01ns. 



